Part A: Time and Frequency domain Analysis

1 Plotting the data and exploring seasonality

We will start of by reading the Gin data and plotting it to observe the general structure of the time series.

setwd("/Volumes/SD/GoogleDrive/0.DA/Assignments/5")
library(TSA)
#Handle timeseries objects in R

## Get the beer data and plot it
gin = read.csv("gin.csv",header=TRUE,dec=".",sep=",")
#convert to a timeseries object: specify the "start" and the frequency of repeating measurements
gin = ts(gin[,2],start=1968,freq=12)
plot(gin,lw=2,col="blue")

We can notice there is strong seasonality and a not much of a trend.

Now we will plot autocorrelations with different lag values.

acf(gin, lag.max = 20, type = "correlation", plot = T, drop.lag.0 = F, main = "Correlogram of the data")

The time series is not that stationary as autocorrelations change over different values of lags (time) so there are obviously some dependencies and we can’t say all instances have the same expected value.

We will also try a partial autocorrelation factor correlogram to inspect the unique contribution of each time-lag.
We will do this because previous values can often be correlated with their previous values so we have to “partial out” those effects. So using partial autocorrelation we can avoid the fact that normal autocorrelation takes into consideration long-term trends.

pacf(gin, lag.max = 20, plot = T, main = "Correlogram of the data")

Again we see that there are significant correlations with values of 12 months ago and 1-2 months ago.
We can hardly say that the time series is not stationary.

We will now plot the STL Decomposition to see how Trend, Seasonality and Remainder look.

gin_stl <- stl(gin, s.window = "periodic")
plot(gin_stl, main = "STL Decomposition of Gin Data")

We can observe there is strong seasonality each year but the trend is pretty weak.

We will check autocorrelation and partial autocorrelation of the remainder (de-trended and de-seasoned data) and compare it with the original model.

res = gin_stl$time.series[,3]
par(mfrow = c(1,2))
acf(gin, main = "Correlogram of gin")
acf(res, main = "Correlogram of remainder")

par(mfrow = c(1,2))
pacf(gin, main = "Correlogram of gin")
pacf(res, main = "Correlogram of remainder")

We can notice that the residual isn’t much better than the original sequence and still has outlier values so it wouldn’t be characterized as strictly stationary.
We also observed some improvements with low lag values that they are lower in the residual compared to the original data.

2 Fitting an ARIMA model

Now we will fit an ARIMA model to our data.
Using the forecast package it will automatically find the best p,q and d values for our ARIMA model based on the AICC metric we mentioned before.
We will then use the model to predict data until 2030 in the future and also plot the related 95% confidence interval:

library(forecast)
fit <- auto.arima(gin,max.p = 25,max.q = 25, max.d = 25,seasonal = TRUE,ic = 'aicc')


plot(gin,xlim=c(1968,2030),lw=2,col="blue")
lines(predict(fit,n.ahead=13*12)$pred,lw=2,col="red")
lines(predict(fit,n.ahead=13*12)$pred+1.96*predict(fit,n.ahead=13*12)$se,lw=2,lty="dotted",col="red")
lines(predict(fit,n.ahead=13*12)$pred-1.96*predict(fit,n.ahead=13*12)$se,lw=2,lty="dotted",col="red")

The models forecast (predictions) until 2030 are in a acceptable (although very wide range).
The middle line containing actual predictions seems to reduce the variance of the values (which will occur) with time. In this case there is still seasonality but it’s magnitude is lower than in the actual time-series.
The upper and lower bounds (confidence interval) reasonably reflect the likely values that will occur based on past trends.

Let’s now explore our model to understand the p,q and d coefficients of our model:

summary(fit)
## Series: gin 
## ARIMA(4,0,2)(2,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          ar1     ar2      ar3     ar4      ma1      ma2    sar1    sar2
##       1.3478  0.3138  -0.8478  0.1617  -0.2753  -0.7033  0.3417  0.4254
## s.e.  0.1361  0.2619   0.1475  0.0466   0.1323   0.1299  0.0400  0.0413
##          mean
##       25.7453
## s.e.   0.0636
## 
## sigma^2 estimated as 0.1659:  log likelihood=-313.72
## AIC=647.43   AICc=647.81   BIC=691.4
## 
## Training set error measures:
##                        ME      RMSE       MAE         MPE     MAPE
## Training set 0.0001668324 0.4042789 0.3157413 -0.02415516 1.228126
##                   MASE         ACF1
## Training set 0.2930961 -0.004206352

The shorthand notation used here is:
ARIMA(p, d, q)(P, D, Q)[S]
where p = non-seasonal Auto-Regressive model order, d = non-seasonal differencing, q = non-seasonal Moving Average order,
P = seasonal Auto-Regressive model order, D = seasonal differencing, Q = seasonal Moving Average order,
and S = time span of repeating seasonal pattern.

Specifically for our model, ARIMA(4, 0, 2) represents the non-seasonal part of our model and means we are describing the response variable by combining a 4th order Auto-Regressive model and a 2nd order Moving Average model.

The second set of coefficient (2,0,0) refers to the seasonal part of the model and means we are describing the response variable using a 2nd order Auto-Regressive model.

The last value of 12 represents the number of periods per season (ie. time span of repeating seasonal pattern). From our knowledge of the model this is a very correct calculation and we will use it later.

3 Power Spectral Density

Now we will perform Power Spectral Density analysis on our time series to calculate which periods are the most influential.
Periodogram prodcued by Power Spectral Density analysis is ideal for identifying periodicity in data and estimating the frequency of the period. Because we have already seen strong seasonality on a yearly basis (12 month period) we expect similarly a period of 12 months (even the ARIMA model showed this).

n = length(gin)
Periodogram = Mod(fft(gin-mean(gin)))^2/n
Freq = (1:n -1)/n
Freq = 1 / Freq
plot(x=Freq[1:70], y=Periodogram[1:70], type='h', lwd=3, ylab="Spectral density",xlab="Period (months)")
title("Gin periodogram")

We see that most dense periods are of less than 100 months so we will filter the plot above to only that size.

Let’s now see more precise values by limiting the range of periods on x-axis:

plot(x=Freq[1:70], y=Periodogram[1:70], type='h', lwd=3, ylab="Spectral density",xlab="Period (months)", xlim=c(1, 100))
title("Gin periodogram")

We can see largest densities around 10-15, 40-45 and 60 months.

Let’s see their exact values:

y = cbind(1:70, Freq[1:70], Periodogram[1:70]);
head(y[order(-y[,3]),2:3],3)
##          [,1]      [,2]
## [1,] 12.00000 249.08016
## [2,] 42.85714  44.59593
## [3,] 60.00000  40.12050

We notice the critical period of 12 months is the most important one with spectral density of 249.
The second most critical frequency is 43 months (42.86) with density of 44.6.
The third most critical frequency is 60 months with a density of 40.12.

So the three critical period are 12, 43 and 60 months.

Part B: Image analysis in frequency field

We will start of by loading the images and turning them into gray scale:

library(EBImage)
library(pracma)
library(ptw)

image <- readImage('image.png')
image2 <- readImage('lincoln.jpg')

grayImage <- EBImage::channel(image,"gray")
grayImage2 <- EBImage::channel(image2,"gray")

par(mfrow = c(1,2))
display(grayImage, method="raster")
display(grayImage2, method="raster")

4 Applying DFT

We will now apply a discrete Fourier transform to our images(s):

# apply Fourier transform
image_dft <- fft(grayImage)

#compute magnitudes
image_magn <- sqrt(abs(image_dft))

# center matrix
center_matrix <- function(m) {
  dims <- dim(m)
  m1 <- cbind(m[,(dims[2]/2+1):dims[2]],m[,1:(dims[2]/2)])
  rbind(m1[(dims[1]/2+1):dims[1],],m1[1:(dims[1]/2),])
}

#shift frequencies to the middle (like fftshift)
image_centered <- center_matrix(image_magn)

# get the log value
image_centered_log = normalize(log(image_centered))

# plot
par(mfrow = c(1,2))
display(image, method="raster")
title("Original photo")
display(image_centered_log, method="raster")
title("DFT")

Because there is horizontal text on the image but which is mainly low frequency we can see that the signal in the middle is strongest.
Since there is also vertical (vertically alligned) text the horizontal lines are strong but also mainly in the middle becuase of low frequency.
There are also strong diagonal lines in the letters and that is also noticeable in the DFT.

We will apply DFT again to a different image:

#Apply Fourier Transform
image_dft <- fft(grayImage2)

#compute magnitudes
image_magn <- sqrt(abs(image_dft))

# center matrix
center_matrix <- function(m) {
  dims <- dim(m)
  m1 <- cbind(m[,(dims[2]/2+1):dims[2]],m[,1:(dims[2]/2)])
  rbind(m1[(dims[1]/2+1):dims[1],],m1[1:(dims[1]/2),])
}

#shift frequencies to the middle (like fftshift)
image_centered <- center_matrix(image_magn)

signedlog10 = function(x) {
  ifelse(abs(x) <= 1, 0, sign(x)*log10(abs(x)))
}

# get the log value
image_centered_log = normalize(log(image_centered))

# plot
par(mfrow = c(1,2))
display(grayImage2, method="raster")
title("Original photo")
display(image_centered_log, method="raster")
title("DFT")

In this picture the horizontal and vertical line in the spectrum are somewhat more pronounced than other lines meaning that there is lots of high and low frequency signals in both vertical and horizontal directions. We can also see that there aren’t as many diagonal or such angled lines.
Also because the frequencies are both low and high are across the main lines horizontal and vertical lines, we can notice that intensites in the DFT are high across the whole range.

5 Filtering

a) Low pass filter

Now we will create a function to apply a low pass filter to our data:

# We will use the makeBrush function (of EBImage) to create these.
low_pass_filter = function(image, threshold){
  sizeim = dim(image)[1]
  size = threshold
  w0 = makeBrush(size = size, shape = 'disc')
  
  fill_in <- 1
  if(mod(sizeim,2) == 0){
    fill_in <- 0
  }

  # This creates a disc-sized element of size = threshold, let's pad the rest with 0
  w1=padzeros(w0,((sizeim-size)/2)-fill_in,side="both")
  w1=t(w1)
  w=padzeros(w1,((sizeim-size)/2) - fill_in,side="both")
  #that gives us a sizeim X sizeim filter, I need more zeros to pad
  w=rbind(w,array(0,sizeim - fill_in - 1))
  w=cbind(w,array(0,sizeim - fill_in))
  # center my element such as that the discoid is in the center
  w=center_matrix(center_matrix(w))

  #Remember that filtering in frequency field means just a multiplication!
  #Notice that we apply the filter on the DC-centered-image
  image_flow <- w * as.Image(center_matrix(fft(image)))

  # place the DC component back to its place before the inverse fourier
  image_low  <- fft(center_matrix(image_flow), inverse=TRUE)
  image_out <- normalize(abs(image_low))
  return (image_out)
}

We will now test it on our images at various thresholds:

par(mfrow = c(1,2))
display(low_pass_filter(grayImage,25), method = "raster")
title("Low pass filter = 25")
display(grayImage, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(low_pass_filter(grayImage,55), method = "raster")
title("Low pass filter = 55")
display(grayImage, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(low_pass_filter(grayImage,85), method = "raster")
title("Low pass filter = 85")
display(grayImage, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(low_pass_filter(grayImage,115), method = "raster")
title("Low pass filter = 115")
display(grayImage, method = "raster")
title(main = "Original")

We can notice that images are very blurry and low level details are lost when we apply a low thresholds to our low pass filter.
So the image loses power with very low thresholds.

par(mfrow = c(1,2))
display(low_pass_filter(grayImage2,25), method = "raster")
title("Low pass filter = 25")
display(grayImage2, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(low_pass_filter(grayImage2,55), method = "raster")
title("Low pass filter = 55")
display(grayImage2, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(low_pass_filter(grayImage2,85), method = "raster")
title("Low pass filter = 85")
display(grayImage2, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(low_pass_filter(grayImage2,115), method = "raster")
title("Low pass filter = 115")
display(grayImage2, method = "raster")
title(main = "Original")

We can notice that the lower threshold we apply the more the details we lose.
On the contrary, higher thresholds applied to the low pass filter lead to more details.
Also we notice the general structure of the image remains, although it looks a bit blurry.

b) High pass filter

Now we will create a function to apply a high pass filter to our data:

# We will use the makeBrush function (of EBImage) to create these.
high_pass_filter = function(image, threshold){
  sizeim = dim(image)[1]
  size = threshold
  w0 = makeBrush(size = size, shape = 'disc')
  
  fill_in <- 1
  if(mod(sizeim,2) == 0){
    fill_in <- 0
  }

  # This creates a disc-sized element of size = threshold, let's pad the rest with 0
  w1=padzeros(w0,((sizeim-size)/2)-fill_in,side="both")
  w1=t(w1)
  w=padzeros(w1,((sizeim-size)/2) - fill_in,side="both")
  #that gives us a sizeim X sizeim filter, I need more zeros to pad
  w=rbind(w,array(0,sizeim - fill_in - 1))
  w=cbind(w,array(0,sizeim - fill_in))
  # center my element such as that the discoid is in the center
  w=center_matrix(center_matrix(w))

  # extra line to apply high pass filter !!
  w <- apply(w, 1:2, function(x) x= abs(x -1))

  #Remember that filtering in frequency field means just a multiplication!
  #Notice that we apply the filter on the DC-centered-image
  image_flow <- w * as.Image(center_matrix(fft(image)))

  # place the DC component back to its place before the inverse fourier
  image_low  <- fft(center_matrix(image_flow), inverse=TRUE)
  image_out <- normalize(abs(image_low))
  return (image_out)
}

We will now test it on our images at various thresholds:

par(mfrow = c(1,2))
display(high_pass_filter(grayImage,25), method = "raster")
title("High pass filter = 25")
display(grayImage, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(high_pass_filter(grayImage,55), method = "raster")
title("High pass filter = 55")
display(grayImage, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(high_pass_filter(grayImage,85), method = "raster")
title("High pass filter = 85")
display(grayImage, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(high_pass_filter(grayImage,115), method = "raster")
title("High pass filter = 115")
display(grayImage, method = "raster")
title(main = "Original")

We can notice that images experience a bit of “ringing” with lower thresholds.
On the other hand with higher values we get a relatively sharp look at the edges of the picture but with some loss of detail (structure that fills the letters).

par(mfrow = c(1,2))
display(high_pass_filter(grayImage2,25), method = "raster")
title("High pass filter = 25")
display(grayImage2, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(high_pass_filter(grayImage2,55), method = "raster")
title("High pass filter = 55")
display(grayImage2, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(high_pass_filter(grayImage2,85), method = "raster")
title("High pass filter = 85")
display(grayImage2, method = "raster")
title(main = "Original")

par(mfrow = c(1,2))
display(high_pass_filter(grayImage2,115), method = "raster")
title("High pass filter = 115")
display(grayImage2, method = "raster")
title(main = "Original")

Applying the high pass filter to Lincoln image leads to more loss main in the main picture structure.
We notice mainly the small details and especially the edges.
The higher the threshold also the fewer edges we notice.

6 Mixing the images

Let’s first load the two new images and see how they look:

image <- readImage('artist1.jpg')
image2 <- readImage('artist2.jpg')

grayImage <- EBImage::channel(image,"gray")
grayImage2 <- EBImage::channel(image2,"gray")

par(mfrow = c(1,2))
display(grayImage, method="raster")
display(grayImage2, method="raster")

I think the first image is more suitable for a low pass filter while the second one is better for a low pass filter.

We will plot the first with a low pass filter

par(mfrow = c(1,2))
display(low_pass_filter(grayImage,105), method = "raster")
title("Low pass filter = 105")
display(grayImage, method = "raster")
title(main = "Original")

We can see that low pass filter is good for this picture as it has solid structure and doesn’t have many edges which would be lost after applying a low pass filter.

We are also going to plot the second image with a high pass filter applied:

par(mfrow = c(1,2))
display(high_pass_filter(grayImage2,65), method = "raster")
title("High pass filter = 65")
display(grayImage2, method = "raster")
title(main = "Original")

We can see that high pass filter is good for this picture as it has many edges which remain after the high pass filter is applied.

Mixing the images

Now we will mix the two images while applying a low pass filter to the first image and a high pass filter to the second image in the frequency domain.
We do this by first transforming and filtering each image, then adding them in the frequency domain and lastly converting their mix back to the spatial domain.
We will create a function for that which takes in inputs images 1 and 2, a threshold for low pass filter applied to image 1 and a threshold for high pass filter applied to image 2:

# We will use the makeBrush function (of EBImage) to create these.
combined_filter = function(image1, low_pass_threshold1, image2, high_pass_threshold2){
  # Process image 1
  sizeim1 = dim(image1)[1]
  size1 = low_pass_threshold1
  w01 = makeBrush(size = size1, shape = 'disc')
  fill_in1 <- 1
  if(mod(sizeim1,2) == 0){
    fill_in1 <- 0
  }
  w1_1=padzeros(w01,((sizeim1-size1)/2) - fill_in1,side="both")
  w1_1=t(w1_1)
  w1=padzeros(w1_1,((sizeim1-size1)/2) - fill_in1,side="both")
  w1=rbind(w1,array(0,sizeim1 - fill_in1 - 1))
  w1=cbind(w1,array(0,sizeim1 - fill_in1))
  w1=center_matrix(center_matrix(w1))

  # low pass filter on image 1
  image1_flow <- w1 * as.Image(center_matrix(fft(image1)))
  
  
  # Process image 2
  sizeim2 = dim(image2)[1]
  size2 = high_pass_threshold2
  w02 = makeBrush(size = size2, shape = 'disc')
  fill_in2 <- 1
  if(mod(sizeim2,2) == 0){
    fill_in2 <- 0
  }
  w1_2=padzeros(w02,((sizeim2-size2)/2) - fill_in2,side="both")
  w1_2=t(w1_2)
  w2=padzeros(w1_2,((sizeim2-size2)/2) - fill_in2,side="both")
  w2=rbind(w2,array(0,sizeim2 - fill_in2 - 1))
  w2=cbind(w2,array(0,sizeim2 - fill_in2))
  w2=center_matrix(center_matrix(w2))
  
  # extra line to apply high pass filter !!
  w2 <- apply(w2, 1:2, function(x) x= abs(x -1))
  # high pass filter on image 2
  image2_fhigh <- w2 * as.Image(center_matrix(fft(image2)))
  
  
  #combine low and high filtered images
  image_combined = image1_flow + image2_fhigh
  
  # place the DC component back to its place before the inverse fourier
  image_combined  <- fft(center_matrix(image_combined), inverse=TRUE)
  image_out <- normalize(abs(image_combined))
  return (image_out)
}

Now let’s see how it works on our images:

par(mfrow = c(1,3))

display(grayImage, method="raster")
title(main = "Image1")

display(grayImage2, method="raster")
title(main = "Image2")

display(combined_filter(grayImage,55,grayImage2,105), method = "raster")
title("LPF(55) on Im1 + HPF(105) on Im2")

We will do this multiple times for different thresholds and then comment on the results.

display(combined_filter(grayImage,45,grayImage2,135), method = "raster")
title("LPF(45) on Image 1 + HPF(135) on Image 2")

display(combined_filter(grayImage,75,grayImage2,75), method = "raster")
title("LPF(75) on Image 1 + HPF(75) on Image 2")

display(combined_filter(grayImage,65,grayImage2,125), method = "raster")
title("LPF(65) on Image 1 + HPF(125) on Image 2")

display(combined_filter(grayImage,105,grayImage2,65), method = "raster")
title("LPF(105) on Image 1 + HPF(65) on Image 2")

The last image maybe even looks the best because we applied a high threshold for the low pass filter so the main structure of it doesn’t look very blurry. We also applied a low threshold for the high pass filter so many details were left from that image which is rich in details and edges.
In total there is both a lot of structure from image1 and details from image2.

Generally, we can notice that:
- in the mixed image we keep the main objects (big continuous things) from image 1 to which low pass filter was applied
- the mixed image has mainly details (edges and such) from the image 2 to which a hi gh pass filter was applied

Mixing in spatial vs frequency domain

Whether we add the two images in spatial or frequency field does not really matter.

You can add two signals together in the spatial domain, apply a Fourier transform and get the same result as if adding the Fourier transforms of the two signals together.

To make sure the results would be the same you’d have to match the frequency threshold to spatial intensities and this is somewhat complex.

Let’s try out how this would (approximately) work and see the results.

We will compare adding the two images in spatial domain vs adding them in frequency domain.
To add them in frequency domain we will re-use the function we just defined.
To add the two images in spatial domain we will re-use the high and low pass filter functions we defined a in task 5 (B2).

par(mfrow = c(2,2))

display(grayImage, method="raster")
title(main = "Image1")

display(grayImage2, method="raster")
title(main = "Image2")

display(low_pass_filter(grayImage,105) + high_pass_filter(grayImage2,65), method = "raster")
title(main = "Spatial LPF(105) on Im1 + HPF(65) on Im2")

display(combined_filter(grayImage,105,grayImage2,65), method = "raster")
title("Frequency LPF(105) on Im1 + HPF(65) on Im2")

And another attempt:

par(mfrow = c(2,1))

display(low_pass_filter(grayImage,45) + high_pass_filter(grayImage2,135), method = "raster")
title(main = "Spatial LPF(45) on Im1 + HPF(135) on Im2")

display(combined_filter(grayImage,45,grayImage2,135), method = "raster")
title("Frequency LPF(45) on Im1 + HPF(135) on Im2")

We can see the images are indeed very similar but unfortunately there are small differences in intensities likely due to reasons we mentioned before with matching of spatial intensities to frequency threshold.